library(knitr)
# purl("h2fec_lifetime_MANOVA.Rmd", output = "h2fec_lifetime_MANOVA.R", documentation = 2)

Setup

rerun <- FALSE

set.seed(224819)

library(MCMCglmm)
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 2.2.1     ✔ purrr   0.2.4
## ✔ tibble  1.4.2     ✔ dplyr   0.7.4
## ✔ tidyr   0.8.0     ✔ stringr 1.3.0
## ✔ readr   1.1.1     ✔ forcats 0.3.0
## ── Conflicts ────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(parallel)
h2fec <- read.table("../../Data/Processed/eggs_per_vial.txt",
                    sep = '\t',
                    dec = ".", header = TRUE,
                    stringsAsFactors = FALSE)

# Standardize egg_total
h2fec$egg_total <- as.numeric(scale(h2fec$egg_total))

h2fec$animal <- seq(1, nrow(h2fec))
h2fec$treat <- as.factor(h2fec$treat)

pedigree <- h2fec[, c("animal", "sireid", "damid")]
names(pedigree) <- c("animal", "sire", "dam")
pedigree$animal <- as.character(pedigree$animal)
pedigree$sire <- as.character(pedigree$sire)
pedigree$dam <- as.character(pedigree$dam)
sires <- data.frame(animal = unique(pedigree$sire),
                    sire = NA, dam = NA, stringsAsFactors = FALSE)
dams <- data.frame(animal = unique(pedigree$dam),
                   sire = NA, dam = NA, stringsAsFactors = FALSE)
pedigree <- bind_rows(sires, dams, pedigree) %>%
  as.data.frame()

genet_corr <- tibble(model = character(),
                     HS_STD = numeric(),
                     LY_STD = numeric(),
                     HS_LY = numeric(),
                     n_eff = numeric())

iter <- 6.5e6
burnin <- 5e4
thinning <- 500
chains <- 12

MANOVA analysis

# 19 hours run time on nivalis (6e5 would be reasonable!)
if (rerun) {
  HS <- h2fec %>% 
    filter(treat == "HS") %>% rename(Eggs_HS = egg_total) %>% 
    as.data.frame()
  LY <- h2fec %>% 
    filter(treat == "LY") %>% rename(Eggs_LY = egg_total) %>% 
    as.data.frame()
  STD <- h2fec %>%
    filter(treat == "STD") %>% rename(Eggs_STD = egg_total) %>% 
    as.data.frame()
  
  h2fec_mrg <- full_join(full_join(HS, LY), STD)
  
  prior1 <- list(R = list(V = diag(3) * 1.002, nu = 1.002),
                 G = list(G1 = list(V = diag(3) * 1.002, nu = 0.002)))
  
  priors <- list(prior1)
  prior_names <- c("Tri: V = diag(3) * 1.002, nu = 0.02")
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(cbind(Eggs_STD, Eggs_HS, Eggs_LY) ~ trait - 1,
               random = ~ us(trait):animal,
               rcov = ~ idh(trait):units,
               data = h2fec_mrg,
               prior = prior,
               pedigree = pedigree,
               family = rep("gaussian", 3),
               nitt = iter,
               burnin = burnin,
               thin = thinning)
    }, mc.cores = chains)
    outfile <- paste0("../../Data/Processed/FEC_lifetime_model_prior", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
      sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
             re[ , "traitEggs_HS:traitEggs_HS.animal"])
    
    # STD vs. LY
    LY_STD <- re[ , "traitEggs_LY:traitEggs_STD.animal"] /
      sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
             re[ , "traitEggs_LY:traitEggs_LY.animal"])
    
    # LY vs. HS
    HS_LY <- re[ , "traitEggs_HS:traitEggs_LY.animal"] /
      sqrt(re[ , "traitEggs_LY:traitEggs_LY.animal"] *
             re[ , "traitEggs_HS:traitEggs_HS.animal"])
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   LY_STD = median(LY_STD),
                                   HS_LY = median(HS_LY),
                                   n_eff = mean(n_eff)))
  }
  
  write_csv(genet_corr, path = "../../Data/Processed/Genetic_Correlations_Fecundity.csv")
}

Analyze model

load("../../Data/Processed/FEC_lifetime_model_prior1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)

plot(fe[, 1, drop = FALSE], ask = FALSE)

plot(fe[, 2, drop = FALSE], ask = FALSE)

plot(fe[, 3, drop = FALSE], ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1, drop = FALSE], ask = FALSE)

plot(re[, 2, drop = FALSE], ask = FALSE)

plot(re[, 3, drop = FALSE], ask = FALSE)

autocorr.diag(re)
##           traitEggs_STD:traitEggs_STD.animal
## Lag 0                           1.0000000000
## Lag 500                         0.1301872307
## Lag 2500                        0.0238423060
## Lag 5000                       -0.0009196147
## Lag 25000                      -0.0033017312
##           traitEggs_HS:traitEggs_STD.animal
## Lag 0                          1.0000000000
## Lag 500                        0.3584259603
## Lag 2500                       0.0762444644
## Lag 5000                       0.0119536087
## Lag 25000                     -0.0009957461
##           traitEggs_LY:traitEggs_STD.animal
## Lag 0                          1.0000000000
## Lag 500                        0.2315114560
## Lag 2500                       0.0455266379
## Lag 5000                       0.0016084792
## Lag 25000                     -0.0007429788
##           traitEggs_STD:traitEggs_HS.animal
## Lag 0                          1.0000000000
## Lag 500                        0.3584259603
## Lag 2500                       0.0762444644
## Lag 5000                       0.0119536087
## Lag 25000                     -0.0009957461
##           traitEggs_HS:traitEggs_HS.animal
## Lag 0                          1.000000000
## Lag 500                        0.336434748
## Lag 2500                       0.057350136
## Lag 5000                       0.010662096
## Lag 25000                      0.001899512
##           traitEggs_LY:traitEggs_HS.animal
## Lag 0                         1.0000000000
## Lag 500                       0.3657226020
## Lag 2500                      0.0764135329
## Lag 5000                      0.0085199420
## Lag 25000                    -0.0008647684
##           traitEggs_STD:traitEggs_LY.animal
## Lag 0                          1.0000000000
## Lag 500                        0.2315114560
## Lag 2500                       0.0455266379
## Lag 5000                       0.0016084792
## Lag 25000                     -0.0007429788
##           traitEggs_HS:traitEggs_LY.animal
## Lag 0                         1.0000000000
## Lag 500                       0.3657226020
## Lag 2500                      0.0764135329
## Lag 5000                      0.0085199420
## Lag 25000                    -0.0008647684
##           traitEggs_LY:traitEggs_LY.animal traitEggs_STD.units
## Lag 0                         1.0000000000         1.000000000
## Lag 500                       0.3665963079         0.122943780
## Lag 2500                      0.0866587765         0.024262491
## Lag 5000                      0.0115061204         0.002461673
## Lag 25000                    -0.0002152714        -0.001991352
##           traitEggs_HS.units traitEggs_LY.units
## Lag 0           1.0000000000        1.000000000
## Lag 500         0.1379153141        0.190259390
## Lag 2500        0.0218513151        0.041620629
## Lag 5000        0.0006775605        0.004527489
## Lag 25000       0.0001317325       -0.001041116
effectiveSize(re)
## traitEggs_STD:traitEggs_STD.animal  traitEggs_HS:traitEggs_STD.animal 
##                           90079.97                           49411.80 
##  traitEggs_LY:traitEggs_STD.animal  traitEggs_STD:traitEggs_HS.animal 
##                           64857.14                           49411.80 
##   traitEggs_HS:traitEggs_HS.animal   traitEggs_LY:traitEggs_HS.animal 
##                           54354.90                           48082.96 
##  traitEggs_STD:traitEggs_LY.animal   traitEggs_HS:traitEggs_LY.animal 
##                           64857.14                           48082.96 
##   traitEggs_LY:traitEggs_LY.animal                traitEggs_STD.units 
##                           47568.07                           93396.54 
##                 traitEggs_HS.units                 traitEggs_LY.units 
##                           90928.34                           71209.90
# Concatenate samples
re <- as.mcmc(as.matrix(re))

head(re)
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1 
## End = 7 
## Thinning interval = 1 
##      traitEggs_STD:traitEggs_STD.animal traitEggs_HS:traitEggs_STD.animal
## [1,]                          0.7873924                         0.2284112
## [2,]                          0.5321631                         0.1727901
## [3,]                          0.2046517                        -0.0436210
## [4,]                          0.1470370                         0.0380646
## [5,]                          0.4637641                         0.2837085
## [6,]                          0.5093006                         0.1375689
## [7,]                          0.3839360                         0.1349881
##      traitEggs_LY:traitEggs_STD.animal traitEggs_STD:traitEggs_HS.animal
## [1,]                        0.20062836                         0.2284112
## [2,]                        0.16648431                         0.1727901
## [3,]                        0.03403463                        -0.0436210
## [4,]                        0.05747749                         0.0380646
## [5,]                        0.28755552                         0.2837085
## [6,]                        0.21914867                         0.1375689
## [7,]                        0.17367460                         0.1349881
##      traitEggs_HS:traitEggs_HS.animal traitEggs_LY:traitEggs_HS.animal
## [1,]                       0.12390648                       0.08419042
## [2,]                       0.30457326                       0.16529953
## [3,]                       0.07812724                       0.01443818
## [4,]                       0.08090601                       0.04531889
## [5,]                       0.18517285                       0.18483029
## [6,]                       0.04638687                       0.06680849
## [7,]                       0.05159149                       0.06498238
##      traitEggs_STD:traitEggs_LY.animal traitEggs_HS:traitEggs_LY.animal
## [1,]                        0.20062836                       0.08419042
## [2,]                        0.16648431                       0.16529953
## [3,]                        0.03403463                       0.01443818
## [4,]                        0.05747749                       0.04531889
## [5,]                        0.28755552                       0.18483029
## [6,]                        0.21914867                       0.06680849
## [7,]                        0.17367460                       0.06498238
##      traitEggs_LY:traitEggs_LY.animal traitEggs_STD.units
## [1,]                       0.06357273           0.1973706
## [2,]                       0.10326474           0.2148104
## [3,]                       0.01353443           0.4054543
## [4,]                       0.04068348           0.4818574
## [5,]                       0.18706363           0.3321149
## [6,]                       0.10110361           0.2497468
## [7,]                       0.08591935           0.3466639
##      traitEggs_HS.units traitEggs_LY.units
## [1,]         0.29195323          0.5152214
## [2,]         0.09147268          0.4236806
## [3,]         0.32326398          0.5218261
## [4,]         0.21465276          0.4229175
## [5,]         0.25146149          0.3877337
## [6,]         0.25396660          0.3365428
## [7,]         0.36814486          0.3673639
# STD vs. HS
plot(re[ , "traitEggs_HS:traitEggs_STD.animal"])

plot(re[ , "traitEggs_STD:traitEggs_STD.animal"])

HS_STD <- re[ , "traitEggs_HS:traitEggs_STD.animal"] /
  sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
         re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_STD)

median(HS_STD)
## [1] 0.7482707
HPDinterval(HS_STD)
##           lower     upper
## var1 -0.4594869 0.9992129
## attr(,"Probability")
## [1] 0.95
# STD vs. LY
LY_STD <- re[ , "traitEggs_LY:traitEggs_STD.animal"] /
  sqrt(re[ , "traitEggs_STD:traitEggs_STD.animal"] *
         re[ , "traitEggs_LY:traitEggs_LY.animal"])
plot(LY_STD)

median(LY_STD)
## [1] 0.9525121
HPDinterval(LY_STD)
##          lower    upper
## var1 0.5776618 0.999608
## attr(,"Probability")
## [1] 0.95
# LY vs. HS
HS_LY <- re[ , "traitEggs_HS:traitEggs_LY.animal"] /
  sqrt(re[ , "traitEggs_LY:traitEggs_LY.animal"] *
         re[ , "traitEggs_HS:traitEggs_HS.animal"])
plot(HS_LY)

median(HS_LY)
## [1] 0.7855926
HPDinterval(HS_LY)
##          lower     upper
## var1 -0.468109 0.9988741
## attr(,"Probability")
## [1] 0.95

Plot for poster

library(tidyverse)
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
B <- data_frame(`HS vs. STD` = as.numeric(HS_STD),
                `LY vs. STD` = as.numeric(LY_STD),
                `HS vs. LY` = as.numeric(HS_LY))

colMeans(B)
## HS vs. STD LY vs. STD  HS vs. LY 
##  0.5762542  0.8900298  0.6011895
B %>% gather(Comparison, value) %>%
  ggplot(aes(value, color = Comparison)) +
  geom_line(stat = "density", size = 2) +
  labs(x = "Genetic Correlation", y = "Density") +
  theme(legend.position = c(0.15, 0.85),
        text = element_text(size = 24),
        legend.text = element_text(size = 18))

ggsave(last_plot(), file = "Genetic_Correlation_Plot.pdf",
       width = 8, height = 6)

Model following Ingleby et al. 2013

Setup data as above

V_nu_to_a_b <- function(V, nu) {
  require(tidyverse)
  require(invgamma)
  
  a <- nu / 2
  b <- (nu * V) / 2
  
  p <- tibble(
    x = seq(0, 10, length = 200),
    y = dinvgamma(x, a, b)) %>% 
    ggplot(aes(x, y)) +
    geom_line()
  print(p)
  
  return(list(alpha = a, beta = b))
}
V_nu_to_a_b(1, 0.002)
## Loading required package: invgamma
## Warning: Removed 1 rows containing missing values (geom_path).

## $alpha
## [1] 0.001
## 
## $beta
## [1] 0.001
if (rerun) {
  
  prior1 <- list(R = list(V = 1, nu = 0.002),
                 G = list(G1 = list(V = diag(3) / 3, nu = 0.002)))
  
  prior2 <- list(R = list(V = 1, nu = 0.002),
                 G = list(G1 = list(V = diag(3) * 0.5 * 0.7, nu = 0.002)))
  
  M <- matrix(0.5 * 0.7, 3, 3)
  diag(M) <- 0.5
  prior3 <- list(R = list(V = 1, nu = 0.002),
                 G = list(G1 = list(V = M, nu = 0.002)))
  
  priors <- list(prior1, prior2, prior3)
  prior_names <- c("Sire: V = diag(3) / 3, nu = 0.002",
                   "Sire: V = diag(3) * 0.5 * 0.7, nu = 0.002",
                   "Sire: V = 0.5 diag, 0.35 off-diag, nu = 0.002")
  
  h2fec$egg_total_std <- scale(h2fec$egg_total)
  
  for (ii in 1:length(priors)) {
    prior <- priors[[ii]]
    
    fm <- mclapply(1:chains, function(i) {
      MCMCglmm(egg_total_std ~ treat,
               random = ~ us(treat):sireid,
               data = h2fec,
               prior = prior,
               family = "gaussian",
               nitt = iter,
               burnin = burnin,
               thin = thinning,
               verbose = TRUE)
    }, mc.cores = chains)
    
    outfile <- paste0("sire_model_prior", ii, ".Rda")
    save(fm, file = outfile)
    
    re <- lapply(fm, function(m) m$VCV)
    re <- do.call(mcmc.list, re)
    re <- as.mcmc(as.matrix(re))
    n_eff <- effectiveSize(re)
    
    # STD vs. HS
    HS_STD <- re[ , "treatHS:treatSTD.sireid"] /
      sqrt(re[ , "treatHS:treatHS.sireid"] *
             re[ , "treatSTD:treatSTD.sireid"])
    median(HS_STD)
    
    # STD vs. LY
    LY_STD <- re[ , "treatSTD:treatLY.sireid"] /
      sqrt(re[ , "treatSTD:treatSTD.sireid"] *
             re[ , "treatLY:treatLY.sireid"])
    median(LY_STD)
    
    # LY vs. HS
    HS_LY <- re[ , "treatHS:treatLY.sireid"] /
      sqrt(re[ , "treatLY:treatLY.sireid"] *
             re[ , "treatHS:treatHS.sireid"])
    median(HS_LY)
    
    genet_corr <- bind_rows(genet_corr,
                            tibble(model = prior_names[[ii]],
                                   HS_STD = median(HS_STD),
                                   LY_STD = median(LY_STD),
                                   HS_LY = median(HS_LY),
                                   n_eff = mean(n_eff)))
    
  }
  write_csv(genet_corr, path = "Genetic_Correlations.csv")
  
  kable(genet_corr, digits = 3)
}

Analyze model

load("sire_model_prior1.Rda")

fe <- lapply(fm, function(m) m$Sol)
fe <- do.call(mcmc.list, fe)
plot(fe, ask = FALSE)

# Extract random effects, convert to mcmc.list
re <- lapply(fm, function(m) m$VCV)
re <- do.call(mcmc.list, re)

plot(re[, 1:3], ask = FALSE)

plot(re[, 4:6], ask = FALSE)

plot(re[, 7:9], ask = FALSE)

plot(re[, 10], ask = FALSE)

effectiveSize(re)
##   treatHS:treatHS.sireid   treatLY:treatHS.sireid  treatSTD:treatHS.sireid 
##                 154800.0                 153923.7                 153390.7 
##   treatHS:treatLY.sireid   treatLY:treatLY.sireid  treatSTD:treatLY.sireid 
##                 153923.7                 153548.4                 155817.7 
##  treatHS:treatSTD.sireid  treatLY:treatSTD.sireid treatSTD:treatSTD.sireid 
##                 153390.7                 155817.7                 155773.0 
##                    units 
##                 155616.7
# Concatenate samples
re <- as.mcmc(as.matrix(re))
colnames(re)
##  [1] "treatHS:treatHS.sireid"   "treatLY:treatHS.sireid"  
##  [3] "treatSTD:treatHS.sireid"  "treatHS:treatLY.sireid"  
##  [5] "treatLY:treatLY.sireid"   "treatSTD:treatLY.sireid" 
##  [7] "treatHS:treatSTD.sireid"  "treatLY:treatSTD.sireid" 
##  [9] "treatSTD:treatSTD.sireid" "units"
# ??? Heritability ????
HS <- re[, "treatHS:treatHS.sireid"] / (re[, "treatHS:treatHS.sireid"] + re[, "units"])
median(HS)
## [1] 0.02586669
LY <- re[, "treatLY:treatLY.sireid"] / (re[, "treatLY:treatLY.sireid"] + re[, "units"])
median(LY)
## [1] 0.07186056
STD <- re[, "treatSTD:treatSTD.sireid"] / (re[, "treatSTD:treatSTD.sireid"] + re[, "units"])
median(STD)
## [1] 0.3158403
# STD vs. HS
HS_STD <- re[ , "treatHS:treatSTD.sireid"] /
  sqrt(re[ , "treatHS:treatHS.sireid"] *
         re[ , "treatSTD:treatSTD.sireid"])
plot(HS_STD)

median(HS_STD)
## [1] 0.6664934
HPDinterval(HS_STD)
##           lower     upper
## var1 -0.8625958 0.9993754
## attr(,"Probability")
## [1] 0.95
# STD vs. LY
LY_STD <- re[ , "treatSTD:treatLY.sireid"] /
  sqrt(re[ , "treatSTD:treatSTD.sireid"] *
         re[ , "treatLY:treatLY.sireid"])
plot(LY_STD)

median(LY_STD)
## [1] 0.932974
HPDinterval(LY_STD)
##          lower     upper
## var1 0.3210507 0.9997023
## attr(,"Probability")
## [1] 0.95
# LY vs. HS
HS_LY <- re[ , "treatHS:treatLY.sireid"] /
  sqrt(re[ , "treatLY:treatLY.sireid"] *
         re[ , "treatHS:treatHS.sireid"])
plot(HS_LY)

median(HS_LY)
## [1] 0.67033
HPDinterval(HS_LY)
##           lower     upper
## var1 -0.8509148 0.9987058
## attr(,"Probability")
## [1] 0.95